﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace RSEngine.Math
{

    /// <summary>
    /// What is Quasi-Monte Carlo?
    /// 
    /// In numerical analysis, a quasi-Monte Carlo method is a method for the computation of
    /// an integral (or some other problem) that is based on low-discrepancy sequences.
    /// 
    /// The Monte Carlo method is a method of taking a series of samples/readings of values
    /// (any kind of values, such as light values, colour values, reflective states) in or
    /// around an area at random, so as to determine the correct actions to take in certain
    /// calculations which usually require multiple sample values to determine overall accuracy,
    /// of those calculations. The Monte Carlo methods tries to be as random as possible, this can
    /// often cause areas that are being sampled to have large irregular gaps in them
    /// (places that are not sampled/read), this in turn can cause problems for certain
    /// calculations (such as shadow calculation).
    /// 
    /// The solution to this was the Quasi-Monte Carlo method.
    /// 
    /// The Quasi-Monte Carlo method is also random, but tries to make sure that the samples/readings
    /// it takes are also better distributed (leaving less irregular gaps in its sample areas) and more
    /// evenly spread across an area. This has the advantage of sometimes leading to more accurate
    /// calculations based on samples/reading. 
    /// </summary>
    public sealed class QMC
    {

        #region Variables

        private const int NUM = 128;
        
        private static readonly int[][] SIGMA = new int[NUM][];
        
        private static readonly int[] PRIMES = new int[NUM];

        #endregion

        #region Constructor

        static QMC()
        {
            {
                // Initializing Faure scrambling tables

                // build table of first primes
                PRIMES[0] = 2;
                for (int i = 1; i < PRIMES.Length; i++)
                    PRIMES[i] = nextPrime(PRIMES[i - 1]);
                int[][] table = new int[PRIMES[PRIMES.Length - 1] + 1][];
                table[2] = new int[2];
                table[2][0] = 0;
                table[2][1] = 1;
                for (int i = 3; i <= PRIMES[PRIMES.Length - 1]; i++)
                {
                    table[i] = new int[i];
                    if ((i & 1) == 0)
                    {
                        int[] prev = table[i >> 1];
                        for (int j = 0; j < prev.Length; j++)
                            table[i][j] = 2 * prev[j];
                        for (int j = 0; j < prev.Length; j++)
                            table[i][prev.Length + j] = 2 * prev[j] + 1;
                    }
                    else
                    {
                        int[] prev = table[i - 1];
                        int med = (i - 1) >> 1;
                        for (int j = 0; j < med; j++)
                            table[i][j] = prev[j] + ((prev[j] >= med) ? 1 : 0);
                        table[i][med] = med;
                        for (int j = 0; j < med; j++)
                            table[i][med + j + 1] = prev[j + med] + ((prev[j + med] >= med) ? 1 : 0);
                    }
                }
                for (int i = 0; i < PRIMES.Length; i++)
                {
                    int p = PRIMES[i];
                    SIGMA[i] = new int[p];
                    Array.Copy(table[p], 0, SIGMA[i], 0, p);
                }
            }
        }

        #endregion

        #region nextPrime

        private static int nextPrime(int p)
        {
            p = p + (p & 1) + 1;
            while (true)
            {
                int div = 3;
                bool isPrime = true;
                while (isPrime && ((div * div) <= p))
                {
                    isPrime = ((p % div) != 0);
                    div += 2;
                }
                if (isPrime)
                    return p;
                p += 2;
            }
        }

        #endregion

        #region riVDC

        public static double riVDC(int bits, int r)
        {
            bits = (bits << 16) | (URShift(bits, 16));
            bits = ((bits & 0x00ff00ff) << 8) | (URShift((bits & unchecked((int)0xff00ff00)), 8));
            bits = ((bits & 0x0f0f0f0f) << 4) | (URShift((bits & unchecked((int)0xf0f0f0f0)), 4));
            bits = ((bits & 0x33333333) << 2) | (URShift((bits & unchecked((int)0xcccccccc)), 2));
            bits = ((bits & 0x55555555) << 1) | (URShift((bits & unchecked((int)0xaaaaaaaa)), 1));
            bits ^= r;
            
            return (double)(bits & unchecked((int)0xFFFFFFFFL)) / (double)0x100000000L;
        }

        #endregion

        #region riS

        public static double riS(int i, int r)
        {
            for (int v = 1 << 31; i != 0; i = URShift(i, 1), v ^= URShift(v, 1))
                if ((i & 1) != 0)
                    r ^= v;

            return (double)r / (double)0x100000000L;
        }

        #endregion

        #region riLP

        public static double riLP(int i, int r)
        {
            for (int v = 1 << 31; i != 0; i = URShift(i, 1), v |= URShift(v, 1))
                if ((i & 1) != 0)
                    r ^= v;
            //UPGRADE_WARNING: Data types in Visual C# might be different.  Verify the accuracy of narrowing conversions. "ms-help://MS.VSCC.v80/dv_commoner/local/redirect.htm?index='!DefaultContextWindowIndex'&keyword='jlca1042'"
            return (double)r / (double)0x100000000L;
        }

        #endregion

        #region halton

        public static double halton(int d, int i)
        {
            // generalized Halton sequence
            switch (d)
            {

                case 0:
                    {
                        i = (i << 16) | (URShift(i, 16));
                        i = ((i & 0x00ff00ff) << 8) | (URShift((i & unchecked((int)0xff00ff00)), 8));
                        i = ((i & 0x0f0f0f0f) << 4) | (URShift((i & unchecked((int)0xf0f0f0f0)), 4));
                        i = ((i & 0x33333333) << 2) | (URShift((i & unchecked((int)0xcccccccc)), 2));
                        i = ((i & 0x55555555) << 1) | (URShift((i & unchecked((int)0xaaaaaaaa)), 1));
                        //UPGRADE_WARNING: Data types in Visual C# might be different.  Verify the accuracy of narrowing conversions. "ms-help://MS.VSCC.v80/dv_commoner/local/redirect.htm?index='!DefaultContextWindowIndex'&keyword='jlca1042'"
                        return (double)(i & unchecked((int)0xFFFFFFFFL)) / (double)0x100000000L;
                    }

                case 1:
                    {
                        double v = 0;
                        double inv = 1.0 / 3;
                        double p;
                        int n;
                        for (p = inv, n = i; n != 0; p *= inv, n /= 3)
                            v += (n % 3) * p;
                        return v;
                    }

                default:
                    break;

            }
            int base_Renamed = PRIMES[d];
            int[] perm = SIGMA[d];
            double v2 = 0;
            double inv2 = 1.0 / base_Renamed;
            double p2;
            int n2;
            for (p2 = inv2, n2 = i; n2 != 0; p2 *= inv2, n2 /= base_Renamed)
                v2 += perm[n2 % base_Renamed] * p2;
            return v2;
        }

        #endregion

        #region mod1

        public static double mod1(double x)
        {
            // assumes x >= 0
            //UPGRADE_WARNING: Data types in Visual C# might be different.  Verify the accuracy of narrowing conversions. "ms-help://MS.VSCC.v80/dv_commoner/local/redirect.htm?index='!DefaultContextWindowIndex'&keyword='jlca1042'"
            return x - (int)x;
        }

        #endregion

        #region generateSigmaTable

        public static int[] generateSigmaTable(int n)
        {
            bool generatedAux = (n & (n - 1)) == 0;
            int[] sigma = new int[n];
            for (int i = 0; i < n; i++)
            {
                int digit = n;
                sigma[i] = 0;
                for (int bits = i; bits != 0; bits >>= 1)
                {
                    digit >>= 1;
                    if ((bits & 1) != 0)
                        sigma[i] += digit;
                }
            }
            return sigma;
        }

        #endregion

        #region URShift

        /// <summary>
        /// Performs an unsigned bitwise right shift with the specified number
        /// </summary>
        /// <param name="number">Number to operate on</param>
        /// <param name="bits">Ammount of bits to shift</param>
        /// <returns>The resulting number from the shift operation</returns>
        public static int URShift(int number, int bits)
        {
            if (number >= 0)
                return number >> bits;
            else
                return (number >> bits) + (2 << ~bits);
        }

        #endregion

    }

}
